﻿module Gauss



open System
open System.IO
open MathNet.Numerics
open MathNet.Numerics.FSharp
open MathNet.Numerics.LinearAlgebra.Double
open System.Diagnostics
open System.Collections.Generic

//struttura dati di appoggio
type GaussianPoint =
  struct 
    val Point   : DenseVector
    val Label   : float
    new (point, label) = { Point = point; Label = label }
  end

//funzioni ausiliarie
let fi (y:float) (points : List<GaussianPoint>) =
  points |> Seq.filter(fun p -> p.Label = y)  
         |> Seq.length 
         |> (fun sum -> sum / points.Count)


let mu (y:float) (points : List<GaussianPoint>) =
  let empty_vector = DenseVector(points.[0].Point.Count)
  let y_points = points |> Seq.filter(fun p -> p.Label = y) 

  y_points |> Seq.fold(fun k p -> k + p.Point) empty_vector
           |> (fun sum -> sum / (y_points |> Seq.length |> float))


let sigma (points : List<GaussianPoint>) (mu_list:List<DenseVector>) =
  let empty_matrix = DenseMatrix(points.[0].Point.Count)
  points |> Seq.mapi(fun i p -> (p,i))
         |> Seq.fold(fun k (p,i) -> k + DenseVector.OuterProduct( (p.Point - mu_list.[int p.Label]), (p.Point - mu_list.[int p.Label]))) empty_matrix
         |> (fun sum -> sum.Divide(points |> Seq.length |> float) :?> DenseMatrix) 

let p_k (input : List<GaussianPoint>) (k : float) =
  let k_points = List(input |> Seq.filter(fun p -> p.Label = k))
  (float k_points.Count) / (float input.Count)


let prob_gauss_ki (observed_point : DenseVector) (k : float) (points : List<GaussianPoint>) (mu_list:List<DenseVector>) (sigma:DenseMatrix) (sigma_inv: LinearAlgebra.Generic.Matrix<float>) (sigma_det:float) =
  
  let scarto = observed_point - mu_list.[int k]
  (1. / (Math.Pow(2.0 * Math.PI, float observed_point.Count / 2.) * Math.Sqrt(sigma_det))) * Math.Exp(- (scarto * sigma_inv * scarto ) / 2.)

let compute_gausssian (input : List<GaussianPoint>) (y_list : seq<float>) =
  let y_list = List(y_list |> Seq.map(fun k -> k, p_k input k))  
  let mu_list = List(seq{ for y in y_list do yield mu (fst y) input})

  let sigma = sigma input mu_list
  let sigma_inv : LinearAlgebra.Generic.Matrix<float> = sigma.Inverse()
  let sigma_det = sigma.Determinant()
  seq{
    for i in [0..149] do
    
      let observed_point : DenseVector = input.[i].Point
      let p_xy =
        List(seq{for (k, p_k) in y_list do yield (prob_gauss_ki observed_point k input mu_list sigma sigma_inv sigma_det) * p_k})
      
      //stampa del risultato dell'osservazione
      let computer_label = 
        if p_xy.[0] > p_xy.[1] && p_xy.[0] > p_xy.[2] then
          Console.WriteLine("0")
          0.
        elif p_xy.[1] > p_xy.[0] && p_xy.[1] > p_xy.[2] then
          Console.WriteLine("1")
          1.
        else
          Console.WriteLine("2")
          2.

      yield (input.[i].Point, input.[i].Label, computer_label)
  }